---

Initial analysis

import numpy as np
import pandas as pd
import math as m
from statsmodels.tsa.stattools import acf
from statsmodels.graphics.tsaplots import plot_acf
from scipy import *
import scipy.linalg


import scaleogram as scg
import matplotlib.pyplot as plt
import warnings

from visuals import *
from my_lib import *
from SSA_lib import SSA

warnings.simplefilter('ignore')

plt.rcParams['figure.figsize'] = 16, 8
plt.rcParams['font.family'] = 'DejaVu Serif'
plt.rcParams['lines.linewidth'] = 2
plt.rcParams['lines.markersize'] = 8
plt.rcParams['xtick.labelsize'] = 16
plt.rcParams['ytick.labelsize'] = 16
plt.rcParams['legend.fontsize'] = 16
plt.rcParams['axes.titlesize'] = 24
plt.rcParams['axes.labelsize'] = 8

def HankelMatrix(X, L):
    
    N = X.shape[0]
    return scipy.linalg.hankel(X[ : N - L + 1], X[N - L : N])

def CartesianToSpherical(point):
    r = np.sqrt(sum(point ** 2))
    n = len(point)
    phi = np.zeros(n - 1)
    
    for i in range(n - 2):
        phi[i] = np.arccos(point[i] / np.sqrt(sum(point[i:] ** 2)))
        
    if point[-1] >= 0:
        phi[n - 2] = np.arccos(point[n - 2] / np.sqrt(point[n - 1] ** 2 + point[n - 2] ** 2))
    else:
        phi[n - 2] = 2 * np.pi - np.arccos(point[n - 2] / np.sqrt(point[n - 1] ** 2 + point[n - 2] ** 2))
        
    return np.hstack((phi, r))


def TrajectoryToSpherical(tr):
    tr_spherical = np.zeros(tr.shape)
    for i, point in enumerate(tr):
        tr_spherical[i] = CartesianToSpherical(point)
    return tr_spherical


def SphericalToCartesian(point):
    phi, r = point[:-1], point[-1]
    n = len(point)
    x = np.zeros(n) 
    cur = r
    
    for i in range(n - 1):
        x[i] = cur * np.cos(phi[i])
        cur *= np.sin(phi[i])

    x[n - 1] = cur
    return x
    

def TrajectoryToCartesian(tr):
    tr_cartesian = np.zeros(tr.shape)
    for i, point in enumerate(tr):
        tr_cartesian[i] = SphericalToCartesian(point)
    return tr_cartesian

def matpower(a, alpha):
    a = (a + a.T) / 2
    eigenValues, eigenVectors = np.linalg.eig(a)
    return np.dot(np.dot(eigenVectors, np.diag(eigenValues ** alpha)), eigenVectors.T)

# Sufficient dimension reduction problem:
# X ~ R^p
# Y ~ R
# Y depends only on beta^T X
# beta: p x r

def directional_regression(x, y, h, r):
    np.random.seed(200)
    # x[i, j] = i'th object, j'th feature

    # number of objects
    n = x.shape[0]
    # dimension to reduce, r is target dimension, r <= p
    p = x.shape[1]

    r_var = np.cov(x.T)

    signrt = matpower(r_var, -0.5)

    #xst = (x - np.mean(x, axis=0)) / np.std(x, axis=0)
    xst = np.dot(x - np.mean(x, axis=0), signrt)

    #bin_count = 10
    hist = np.histogram(y, bins=h)
    y_bins = np.fmin(np.digitize(y, hist[1]), h) - 1

    prob = np.zeros(h)
    for i in range(h):
        prob[i] = np.sum(y_bins == i) / n
        
    vxy = np.zeros((p, p, h))
    exy = np.zeros((h, p))
    for i in range(h):
        vxy[:,:,i] = np.cov(xst[y_bins == i, :].T)
        exy[i,:] = np.mean(xst[y_bins == i, :], axis=0)

    mat1 = np.zeros((p, p))
    mat2 = np.zeros((p, p))

    for i in range(h):
#         print("i =", i)
        A = exy[i,:][:,np.newaxis].dot(exy[i,:][:,np.newaxis].T)
        B = vxy[:,:,i] + A
        
#         print('vxy:', vxy[:,:,i])
#         print("B:", B)
        mat1 += prob[i] * B.dot(B)
        mat2 += prob[i] * A

    out = 2 * mat1 + 2 * mat2.dot(mat2) + 2 * np.sum(np.diag(mat2))*mat2 - 2 * np.eye(p)
    
    
#     print("prob:", prob)
#     print("mat1:", mat1)
    
    eigenValues, eigenVectors = np.linalg.eig(out)

    idx = eigenValues.argsort()[::-1]
    eigenValues = eigenValues[idx]
    eigenVectors = eigenVectors[:,idx]

    output = np.zeros((n, r))

    for i in range(r):
        u = eigenVectors[i,:]
#         print('Eigen vector {}'.format(u))
        output[:, i] = u.dot(x.T)

    return output

def ApplyDR(trajectory, r=2, h=30):
    y = trajectory[:, -1]
    X = trajectory[:, :-1]
    
    pred = directional_regression(X, y, h, r)
    return pred

def add_free(X):
    return np.hstack([X, np.ones(shape=(X.shape[0], 1))])

def MakeStationary(ts):
    new_ts = []
    for i in range(1, len(ts)):
        new_ts.append(ts[i] - ts[i - 1])
    return np.array(new_ts)

def sph_harm_to_shtools(c):
    """ Convert our list format for the sph harm coefficients/harmonics to pyshtools (2, n, n) format. """
    n = len(c)
    real = is_real_sft(c)
    dim1 = 1 if real else 2
    out = np.zeros((dim1, n, n, *c[0][0].shape)) + 0j
    for l, cc in enumerate(c):
        cc = np.array(cc)
        if not real:
            m_minus = cc[:l][::-1]
            m_plus = cc[l:]
        else:
            m_minus = np.array([])
            m_plus = cc

        # we get warnings here when using reals
        if m_minus.size > 0:
            out[1, l, 1:l+1, ...] = m_minus
        out[0, l, :l+1, ...] = m_plus

    return out

def sph_sample(n, mode='DH'):
    """ Sample grid on a sphere.
    Args:
        n (int): dimension is n x n
        mode (str): sampling mode; DH or GLQ
    Returns:
        theta, phi (1D arrays): polar and azimuthal angles
    """
    assert n % 2 == 0
    j = np.arange(0, n)
    if mode == 'DH':
        return j*np.pi/n, j*2*np.pi/n
    elif mode == 'ours':
        return (2*j+1)*np.pi/2/n, j*2*np.pi/n
    elif mode == 'GLQ':
        from pyshtools.shtools import GLQGridCoord
        phi, theta = GLQGridCoord(n-1)
        # convert latitude to [0, np.pi/2]
        return np.radians(phi+90), np.radians(theta)
    elif mode == 'naive':
        # repeat first and last points; useful for plotting
        return np.linspace(0, np.pi, n), np.linspace(0, 2*np.pi, n)

def DHaj(n, mode='DH'):
    """ Sampling weights. """
    # Driscoll and Healy sampling weights (on the phi dimension)
    # note: weights depend on the chosen grid, given by sph_sample
    if mode == 'DH':
        gridfun = lambda j: np.pi*j/n
    elif mode == 'ours':
        gridfun = lambda j: np.pi*(2*j+1)/2/n
    else:
        raise NotImplementedError()

    l = np.arange(0, n/2)
    a = [(2*np.sqrt(2)/n *
          np.sin(gridfun(j)) *
          (1/(2*l+1) * np.sin((2*l+1)*gridfun(j))).sum())
         for j in range(n)]

    return a


def sph_harm_all(n, as_tfvar=False, real=False):
    """ Compute spherical harmonics for an n x n input (degree up to n // 2)
    Args:
        n (int): input dimensions; order will be n // 2
        as_tfvar (bool): if True, return as list of tensorflow Variables.
        real (bool): if True, return real harmonics
    """
    harmonics = []

    for l in range(n // 2):
        if real:
            minl = 0
        else:
            minl = -l
        row = []
        for m in range(minl, l+1):
            row.append(sph_harm_lm(l, m, n))
        harmonics.append(row)

    else:
        return harmonics
    
def sph_harm_lm(l, m, n):
    """ Wrapper around scipy.special.sph_harm. Return spherical harmonic of degree l and order m. """
    phi, theta = sph_sample(n)
    phi, theta = np.meshgrid(phi, theta)
    f = sph_harm(m, l, theta, phi)

    return f

def is_real_sft(h_or_c):
    """ Detect if list of lists of harmonics or coefficients assumes real inputs (m>0) """
    d = len(h_or_c[1])

    isreal = True if d == 2 else False

    return isreal

---

Choose first 10 steps

Common walking, phone accelerometer

Frequency 0.91$\pm$1 s was found empirically

dt = 459*20
#data = pd.read_csv('data/long_walk_100_acc.csv', delimiter =';', decimal=',')[3815+6*455:3815+6*455+dt]
data = pd.read_csv('data/long_walk_100_acc.csv', delimiter =';', decimal=',')[7009:7009+dt]
#data = pd.read_csv('data/long_walk_100_acc.csv', delimiter =';', decimal=',')[3785:8000]
#data = pd.read_csv('data/home_lin_10_lac.csv', delimiter =';', decimal=',')[:]
frecuency = len(data)/(data['time'].values[-1]-data['time'].values[0])
assert 480 < frecuency < 520

x_acc = ( (data['X_value'].values)**2 + (data['Y_value'].values)**2 + (data['Z_value'].values)**2)**.5
_m = np.mean(x_acc)
x_acc -= _m
t = (data['time'].values).astype(float).reshape([-1,])
t = np.linspace(0,t[-1]-t[0],len(x_acc))

fig = go.Figure()
fig.add_scatter(y = x_acc, mode='lines', name='Sum squares')
fig.show()

t_real = (data['time'].values).astype(float).reshape([-1,])
t_real = t_real - t_real[0]

plt.plot(t_real-t,'.')
plt.grid(which='major',
    color = 'gray', 
    linewidth = 0.8)
plt.minorticks_on()
plt.grid(which='minor', 
    color = 'gray', 
    linestyle = '--',
    linewidth = 0.5)
plt.show()

_m

10.792185704435663

---

Velocity figure

def moving_average(a, n=3):
    ret = np.cumsum(a, dtype=float)
    ret[n:] = ret[n:] - ret[:-n]
    return ret[n - 1:] / n

_t,_x = integ(t,x_acc)

mean_w = 1500

fig = go.Figure()
fig.add_scatter(x = _t[mean_w:], y = (_x)[mean_w:],
                mode='lines', 
                name='Velocity')
fig.add_scatter(x = _t[mean_w:], y = moving_average(_x,mean_w+1),
                mode='lines', 
                name='Velocity')
fig['layout']['yaxis1'].update(title='Velocity[m/s]')
fig['layout']['xaxis1'].update(title='Time[s]')
fig.show()

---

Autocorrelation Function

acf_x = acf(
        x_acc,
        nlags=3000,
        alpha=0.5,
)

fig = go.Figure()
fig.add_scatter(x = np.arange(len(acf_x[0])), y = acf_x[0], mode='lines', name='ACF')
fig.show()

np.arange(0,4000,461)

array([   0,  461,  922, 1383, 1844, 2305, 2766, 3227, 3688])

T_one_step = 460/500
T_one_step

0.92

ax = plot_acf(x_acc,
              lags=len(x_acc)-1,
              use_vlines = False,
              alpha=0.1
         )
plt.tick_params(axis='both', which='major', labelsize=25,length=8, width=4)
plt.grid(which='major',
    color = 'gray', 
    linewidth = 0.8)
plt.minorticks_on()
plt.grid(which='minor', 
    color = 'gray', 
    linestyle = '--',
    linewidth = 0.5)
plt.show()

Fast Fourier transforms¶

Fs = 500

N = len(x_acc)

__Abs = 2.0/N * np.abs(fft(x_acc)[0:N//2])

__FS = np.linspace(0.0, 1.0/2.0*Fs, N//2)

# plotting the spectrum
fig = go.Figure()
fig.add_scatter(x = __FS[0:100],
                y = __Abs[0:100],
                mode='lines',
                name='Magnitude'
               )

temp_fft = fft(x_acc)[0:N//2]

__Ph = np.arctan2(np.imag(temp_fft),np.real(temp_fft))

fig.add_scatter(x = __FS[0:100],
                y = __Ph[0:100],
                mode='markers',
                name='Phase'
               )

fig.add_trace(go.Bar(x= __FS[0:100],
                      y= __Ph[0:100],
                        #mode='markers',
                        name='Phase'
                    )
             )


fig['layout']['yaxis1'].update(title='Magnitude')
fig['layout']['xaxis1'].update(title='Frec[1/s]')
fig.show()

scales = np.linspace(130,450,500)
ax = scg.cws(t[::1], x_acc, 
             figsize = (12,6),
             scales = scales,
             wavelet = 'cmor6-1.5',
             coi=False
            )
ax.set_ylabel('Scale [s]',fontsize = 15)
ax.set_xlabel('Time [s]',fontsize = 15)
plt.show()

Application of SSA

accel_ssa = SSA(x_acc, 500)

plt.figure(figsize=(11, 11))
accel_ssa.plot_wcorr(max=30)
plt.title("W-Correlation for Walking Time Series");

n_fig = 1000
fig = go.Figure()

fig.add_scatter(x = t[:n_fig], y = x_acc[:n_fig], mode='lines', name='Real')

fig.add_scatter(x = t[:n_fig], y = accel_ssa.reconstruct(slice(0,5))[:n_fig], mode='lines', name='SSA')


x = (4.6*np.cos(2*np.pi*2.18*(t[:n_fig])-.24) + 2.39*np.cos(2*np.pi*3.27*(t[:n_fig])-2.46))

fig.add_scatter(x = t[:n_fig], y = x, mode='lines', name='Main frec')

fig.show()

Amplitude modulation

n_fig = 1000
fig = go.Figure()

x = (4.6*np.cos(2*np.pi*2.18*(t[:n_fig])-.24) + 2.39*np.cos(2*np.pi*3.27*(t[:n_fig])-2.46))/np.cos(2*np.pi*2.18*(t[:n_fig])-.24)

fig.add_scatter(x = t[:n_fig], y = x, mode='lines', name='Main frec')

fig.show()

Fs = 500

N = len(x)

__Abs = 2.0/N * np.abs(fft(x)[0:N//2])

__FS = np.linspace(0.0, 1.0/2.0*Fs, N//2)

# plotting the spectrum
fig = go.Figure()
fig.add_scatter(x = __FS[0:30],
                y = __Abs[0:30],
                mode='lines',
                name='Magnitude'
               )

temp_fft = fft(x)[0:N//2]

__Ph = np.arctan2(np.imag(temp_fft),np.real(temp_fft))

fig.add_scatter(x = __FS[0:30],
                y = __Ph[0:30],
                mode='markers',
                name='Phase'
               )

fig.add_trace(go.Bar(x= __FS[0:30],
                     y= __Ph[0:30],
                        #mode='markers',
                        name='Phase'
                    )
             )


fig['layout']['yaxis1'].update(title='Magnitude')
fig['layout']['xaxis1'].update(title='Frec[1/s]')
fig.show()

---

fig = go.Figure()
fig.add_scatter(x = t[:], y = accel_ssa.reconstruct(slice(0,4))[:], mode='lines', name='slice')
fig.show()

import scaleogram as scg
scales = np.linspace(130,450,500)
ax = scg.cws(t[::1], accel_ssa.reconstruct(slice(0,4))[::1], 
             figsize = (12,6),
             scales = scales,
             wavelet = 'cmor3-1.5',
             coi=False
            )
ax.set_ylabel('Scale [s]',fontsize = 15)
ax.set_xlabel('Time [s]',fontsize = 15)
plt.show()

x_acc_clear = accel_ssa.reconstruct(slice(0,5))

track1, basis1 = phase_track(np.array([x_acc_clear]).T, 600, 2)

plot_phase_track(track1[:])

Explained variation for 2 principal components: [0.40011295 0.3676284 ]
Cumulative explained variationfor 2 principal components: 0.7677413524588336

track2, basis2 = phase_track(np.array([x_acc_clear]).T, 300, 3)

Explained variation for 3 principal components: [0.43827307 0.3808992  0.09921792]
Cumulative explained variationfor 3 principal components: 0.9183901925645711

plot_phase_track(track2)

Sphere_projection(track2, n = 1)

eigenvalues_plot(x_acc,600,1,num_l = 50,MinMax = False)

No handles with labels found to put in legend.

eigenvalues_plot(x_acc_clear,600,1,num_l = 50)

No handles with labels found to put in legend.

track5, basis5 = phase_track(np.array([x_acc_clear]).T, 600, 4)

Explained variation for 4 principal components: [0.40011295 0.3676284  0.11352861 0.11324896]
Cumulative explained variationfor 4 principal components: 0.9945189233866947

---

Variance of timeseries

# нормировка
n = 1
_x = track2[:,0] - track2[:,0].mean()
_y = track2[:,1] - track2[:,1].mean()
_z = track2[:,2] - track2[:,2].mean()

track_norm = np.array([_x,_y,_z]).T


for i in range(len(track_norm)):
    track_norm[i] = track_norm[i]/((track_norm[i,:]**2).sum())**.5

# данные для сферы
_t, p = np.mgrid[0:np.pi:100j, 0:2*np.pi:100j]

x_s = 0.9*np.sin(_t)*np.cos(p)
y_s = 0.9*np.sin(_t)*np.sin(p)
z_s = 0.9*np.cos(_t)

#вычисляем среднюю
n_m = 15
track_mean = track_norm[:461,:]/n_m

for i in range(1, n_m):
    track_mean += track_norm[461*i:461*(i+1),:]/n_m

---

To spherical coordinates

def to_spherical(track_norm):
    r = np.ones(len(track_norm))
    phi = np.arctan2(track_norm[:,1],track_norm[:,0])+np.pi
    theta = np.arccos(track_norm[:,2])
    return r,phi,theta

fig = go.Figure()
r,phi,theta = to_spherical(track_norm)
fig.add_scatter(y = theta,
                x = phi,
                mode='markers',
               )

r_m,phi_m,theta_m = to_spherical(track_mean)
fig.add_scatter(y = theta_m,
                x = phi_m,
                mode='markers',
               )
#for i in range(10):
fig['layout']['yaxis1'].update(title='theta [0,pi]')
fig['layout']['xaxis1'].update(title='phi [0,2pi]')

fig.show()

#вычисляем дисперсию
n_m = 15
theta_VAR = (theta[:461]-theta_m[:461])**2/(n_m-1)

for i in range(1, n_m):
    theta_VAR += (theta[461*i:461*(i+1)] - theta_m[:461])**2/(n_m-1)
    
theta_STD = theta_VAR**.5

theta_STD_UP   = (theta_m + 2*theta_STD)
theta_STD_DOWN = (theta_m - 2*theta_STD)

fig = go.Figure()

fig.add_scatter(y = theta,
                x = phi,
                mode='markers',
               )

fig.add_scatter(y = theta_m,
                x = phi_m,
                mode='markers',
               )

fig.add_scatter(y = theta_STD_UP,
                x = phi_m,
                mode='markers',
               )

fig.add_scatter(y = theta_STD_DOWN,
                x = phi_m,
                mode='markers',
               )



fig['layout']['yaxis1'].update(title='theta [0,pi]')
fig['layout']['xaxis1'].update(title='phi [0,2pi]')

fig.show()

import pyshtools as pysh

from scipy.special import sph_harm

def sph_harm_transform(f, mode='DH', harmonics=None):
    """ Project spherical function into the spherical harmonics basis. """
    assert f.shape[0] == f.shape[1]
    
    sumfun = np.sum
    conjfun = np.conj
    n = f.shape[0]
    
    assert np.log2(n).is_integer()

    if harmonics is None:
        harmonics = sph_harm_all(n)

    a = DHaj(n, mode)

    f = f*np.array(a)[np.newaxis, :]

    real = is_real_sft(harmonics)

    coeffs = []
    for l in range(n // 2):
        row = []
        minl = 0 if real else -l
        for m in range(minl, l+1):
            # WARNING: results are off by this factor, when using driscoll1994computing formulas
            factor = 2*np.sqrt(np.pi)
            row.append(sumfun(factor * np.sqrt(2*np.pi)/n * f * conjfun(harmonics[l][m-minl])))
        coeffs.append(row)

    return coeffs

len(theta_m)

461

f = scipy.interpolate.interp1d(np.linspace(0,2,461),theta_m[:]/np.pi)

int(f(2/64)*64)

32

plt.plot(np.linspace(0,1,64),f(np.linspace(0,1,64)))

[<matplotlib.lines.Line2D at 0x1920611d10>]

n_pic = 64*2
a = np.zeros((n_pic,n_pic))
for j in range(n_pic):
    i = int(f(j/n_pic)*n_pic)
    a[i][j] = 1000
    i = int(f(1+j/n_pic)*n_pic)
    a[i][j] = 1000

a = np.diag(np.ones(32))
a = np.concatenate((a,np.zeros((32,32))),axis = 1)
a = np.concatenate((a,np.zeros((32,64))),axis = 0)

a[:10,:10]

array([[1., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 1., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 1., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 1., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 1., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 1., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 1., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 1., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 1., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 1.]])

a = np.array([[1., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.], [1., 1., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.], [0., 1., 1., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.], [0., 0., 1., 1., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.], [0., 0., 0., 1., 1., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.], [0., 0., 0., 0., 1., 1., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0.], [0., 0., 0., 0., 0., 1., 1., 1., 0., 0., 0., 0., 0., 0., 0., 0.], [0., 0., 0., 0., 0., 0., 1., 1., 0., 0., 0., 0., 0., 0., 0., 0.], [0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.], [0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.], [0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.], [0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.], [0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.], [0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.], [0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.], [0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]])

a = np.zeros((128,128))
a[:16,16]

array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.])

a[:32,16] = np.ones(32)

plt.imshow(a)

<matplotlib.image.AxesImage at 0x1926aea990>

a = np.zeros((64,64)) a[20] = 100

import functools
import numpy as np
from scipy.special import sph_harm

w = sph_harm_to_shtools(sph_harm_transform(a.T))
x = pysh.SHCoeffs.from_array(real(w))

grid = x.expand()
fig, ax = grid.plot(show=True)

Directional regression

x = np.array([x_acc]).T #np.array([x_acc_clear]).T
mean = x.mean()
std = x.std()
x = (x - mean) / std

track_dr, basis_dr = phase_track(x, 600, 3)

Explained variation for 3 principal components: [0.27456685 0.25186656 0.07912414]
Cumulative explained variationfor 3 principal components: 0.6055575499884778

track_dr_4, basis_dr_4 = phase_track(x, 600, 100)

Explained variation for 100 principal components: [2.74566846e-01 2.51866562e-01 7.91241424e-02 7.83881178e-02
 2.06826673e-02 2.05908599e-02 1.78290698e-02 1.77259356e-02
 1.64920331e-02 1.57026035e-02 1.34263741e-02 1.28837669e-02
 1.22556213e-02 1.22065841e-02 1.08327934e-02 1.06099542e-02
 9.20022773e-03 9.16530095e-03 8.79397114e-03 8.78186844e-03
 8.17032069e-03 8.12591134e-03 7.78016367e-03 7.74974905e-03
 5.95943992e-03 5.87073713e-03 4.40023254e-03 4.38051019e-03
 3.60113535e-03 3.56157514e-03 2.74327079e-03 2.74022737e-03
 2.28908281e-03 2.20516614e-03 2.18840771e-03 2.17081049e-03
 1.69559336e-03 1.68434584e-03 1.66349289e-03 1.66004707e-03
 1.48197974e-03 1.39607735e-03 1.24642537e-03 1.18407570e-03
 1.04735109e-03 9.68826652e-04 9.49223287e-04 9.02499851e-04
 8.33442588e-04 7.58516276e-04 6.94908453e-04 6.28223645e-04
 5.67841975e-04 4.52789143e-04 3.45668674e-04 3.43015476e-04
 3.29918033e-04 2.56110592e-04 2.17160328e-04 1.71220466e-04
 1.64162321e-04 1.52119803e-04 1.41654329e-04 1.33110476e-04
 1.01451232e-04 8.31894882e-05 8.05822753e-05 7.80043274e-05
 7.77737844e-05 7.04398612e-05 5.75232553e-05 5.68198690e-05
 5.68066015e-05 5.49095855e-05 5.30157567e-05 5.16506155e-05
 4.78294303e-05 4.51232609e-05 4.32132880e-05 4.32120115e-05
 4.25527562e-05 4.24881775e-05 4.13646419e-05 3.96175251e-05
 3.87047812e-05 3.80925492e-05 3.78144976e-05 3.31164415e-05
 3.30791020e-05 2.91982477e-05 2.55768356e-05 2.45914713e-05
 2.32655624e-05 2.06927615e-05 1.99937848e-05 1.90529793e-05
 1.56518071e-05 1.40033358e-05 1.31106698e-05 1.21688398e-05]
Cumulative explained variationfor 100 principal components: 0.9986975216681793

X = HankelMatrix(x_acc_clear,300)

pred_X = ApplyDR(TrajectoryToSpherical(track_dr), r = 2)

spherical_tr = TrajectoryToSpherical(track_dr)

w_dr = np.linalg.lstsq(add_free(pred_X), spherical_tr[:, -1])
y_pred_dr = add_free(pred_X).dot(w_dr[0])

dr_trajectory = spherical_tr
dr_trajectory[:, -1] = y_pred_dr
dr_trajectory = TrajectoryToCartesian(dr_trajectory)

fig_2 = go.Figure()


fig_2.add_trace(go.Scatter3d(x=track_dr[:,0][::n],
                             y=track_dr[:,1][::n],
                             z=track_dr[:,2][::n],
                             #mode='markers',
                             marker=dict(
                                         size=2,
                                         line=dict(
                                                    width=0.1
                                                  )

                                         ),
                             name='trajectory'
                            )
                )
fig_2.layout.template = 'plotly_white'
fig_2.show()

fig_2 = go.Figure()
fig_2.add_trace(go.Scatter3d(x=dr_trajectory[:,0][::n],
                             y=dr_trajectory[:,1][::n],
                             z=dr_trajectory[:,-1][::n],
                             #mode='markers',
                             marker=dict(
                                         size=2,
                                         line=dict(
                                                    width=0.1
                                                  )

                                         ),
                             name='dr_trajectory'
                            )
                )


"""fig_2.add_trace(go.Scatter3d(x=[track_dr[0, 0]],
                           y=[track_dr[0, 1]],
                           z=[track_dr[0, 2]],
                           mode='markers',
                           marker_size=10,
                           marker_color='rgba(255, 10, 0, .7)',
                           name='Start point'))

fig_2.add_trace(go.Scatter3d(x=[track_dr[-1, 0]],
                           y=[track_dr[-1, 1]],
                           z=[track_dr[-1, 2]],
                           mode='markers',
                           marker_size=10,
                           marker_color='rgba(10, 250, 250, .7)',
                           name='End point'))"""

fig_2.layout.template = 'plotly_white'
fig_2.show()

Sphere_projection(dr_trajectory[:,[0,150,-1]], n = 1)

Sphere_projection(track_dr, n = 1)